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Abstract 

We describe in detail a simple and efficient Green Function Monte Carlo 
technique for computing both the ground state energy and the ground state 
properties by the "forward walking" scheme. The simplicity of our reconfig- 
uration process, used to maintain the walker population constant, allows to 
control any source of systematic error in a rigorous and systematic way. We 
apply this method to the Heisenberg model and obtain accurate and reliable 
estimates of the ground state energy, the order parameter and the static spin 
structure factor S(q) for several momenta. For the latter quantity we also 
find very good agreement with available experimental data on the La2CuO^ 

antifer r omagnet . 
75.10.Jm,75.40.Mg,75.30.Ds 
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I. INTRODUCTION 



After almost one decade from the discovery of High-Tc superconductivity we have cer- 
tainly understood much more about magnetism rather than superconductivity. In particular 
since almost all the stechiometric compounds of High-Tc superconductors are good antifer- 
romagnets, well described by the two dimensional Heisenberg model (HM), from the very 
beginning a strong numerical effort has been devoted to the simulation of this model. 
The HM is defined by the following Hamiltonian: 

h j = jy: SiSj (i) 

<i,j> 

where the spin one half vectors Si satisfy S 2 = 3/4 and J is the nearest neighbor antiferro- 
magnetic superexchange coupling, connecting nearest neighbor pairs < i,j >. Henceforth 
periodic boundary conditions are assumed in a finite square lattice with N a = I x I sites. 

Although a rigorous proof that this model has long range antiferromagnetic order in two 
spatial dimension is still lacking, there is a general consensus that long range order exists 
even in this interesting case. In other words its properties should be very well understood 
by the simple spin-wave theory, which assumes long range antiferromagnetic order in the 
ground state. || 

In this work we give accurate ground state properties of the HM using a new and more 
efficient version of the Green Function Monte Carlo (GFMC) technique, first applied on a 
lattice by Trivedi and Ceperley, some years ago @. 

With the present scheme we also estimate the ground state energy of the HM on the 
square lattice to be —0.669442 J ± 0.000026J slightly different, but more accurate of the 
previous GFMC estimates. Analogously we obtain for the antiferromagnetic order parameter 
the value m = 0.3077 ± 0.0004, consistent with other numerical estimates but with a very 
accurate control of the finite size effects. We discuss also our results for the static spin 
structure factor in view of the recently proposed theory for the finite size scaling in a quantum 
antiferromagnet. Q In particular we verify that, as a remarkable prediction of the theory, 
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the small q behavior of this function behaves as S(q) = X C M> where x is the magnetic spin 
susceptibility and c the spin wave velocity of the HM. This relation is particularly important 
as this function is experimentally detectable in neutron scattering experiments. || 

Let us discuss now the technical part of our work, which is based on the GFMC, as we 
have mentioned before. As is well known this technique allows to sample statistically the 
ground state of a many body Hamiltonian if by a set of walkers (wi,Xi) which represent 
vectors of a large (or even infinite) Hilbert space. The set of all configurations x spans 
a normalized and complete basis. The aim of this approach is to sample statistically the 
ground state of H, by a large population of walkers. || As it will be described later on, 
in the finite dimensional case, the GFMC on a lattice is based on a statistical application 
of the Hamiltonian matrix-vector product w^x^ — > (-H)wiXi to the walker configurations 
{wx} i , thus filtering out, after many iterations the desired population distribution for 
the ground state. In this statistical iteration however, the walker weights Wi increase or 
decrease exponentially so that after a few iterations most of the walkers have negligible 
weights and some kind of reconfiguration becomes necessary to avoid large statistical errors. 
The process to eliminate the irrelevant walkers or generate copies of the important ones 
is called "branching". This scheme is in principle exact only if the population of walkers 
is let increase or decrease without any limitation. Any reconfiguration of the population 
size may in fact introduce some spurious correlation between the walkers that may affect 
the statistical sampling of the ground state. In practice for long simulation it is always 
necessary to control the population size, as, otherwise, one easily exceeds the maximum 
allowed computer memory. This control of the walker population size may introduce some 
kind of bias that vanishes quite slowly for infinite number of walkers. |6|]7]] In this case only 
by performing several runs with different numbers of walkers one may in principle estimate 
the size of the residual bias. 



Following the Hetherington's work ||11[| , we define here an efficient reconfiguration process 
at fixed number M of walkers, with a rigorous control of the bias and without need of the 
conventional branching scheme. 



In the last sections we present the results obtained for the HM up to N a = 16 x 16 
lattice size, together with some numerical test on a small N a = 4 x 4 lattice where an 
accurate numerical solution is available by exact diagonalization. Previous calculations on 
this model, using the Green Function Monte Carlo technique, were performed either without 
correcting the bias and controlling it for small lattices with a large number of walkers, 
or by correcting the bias in a way, which is probably correct, but it is not possible to prove 
rigorously. || 



II. THE GFMC TECHNIQUE 

In the following sections we describe in detail how to evaluate the maximum eigenvector 
of a matrix H^i x with all positive definite matrix elements, using a stochastic approach. 
Clearly in any physical problem, described by an Hamiltonian H, the most interesting 
eigenvalue is the lowest one: the ground state energy. This is however just a matter of 
notations, as the ground state of H represents the maximum excited state of —H. In the 
following, for simplicity, we assume to change the sign of H so that the physical ground 
state is, in this notations, the maximum eigenvector of H . More important is instead the 
restriction of positive definite matrix elements H^j x , which drastically constrains the class 
of Hamiltonians that can be treated with this method, without facing the old but still 
unsolved "sign problem". Whenever the Hamiltonian has matrix elements with arbitrary 
sign schemes like "Fixed nodes approximation", and their recent developments to finite 
lattices, are possible within the GFMC method. PJlO| Of course if negative signs occur only 



in the diagonal elements of if a simple change of the hamiltonian H — > H x i x + \8 x t x , will 
not change the ground state but the Hamiltonian will satisfy the condition H x ' tX > for 
large enough shift A. For instance the Heisenberg hamiltonian can be easily casted in the 
previous form as it had been previously shown. || 

From a general point of view the ground state of H can be obtained by applying the well 
known power method: 
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|^o >= lim H L \il) T > 



L — > oo 



(2) 



where the equality holds up to (infinite) normalization, and \ipr > is a trial state non 
orthogonal to the ground state \ipQ>- 

In the following a simple stochastic approach is described for evaluating the state 
H M \tpT >■ To this purpose we define a basic element of this stochastic approach : the 
so called walker. A walker is determined by an index x corresponding to a given element 
\x > of the chosen basis and a weight w. The walker "walks" in the Hilbert space of the 
matrix H and assumes a configuration wx according to a given probability distribution 
P(w, x). 

The task of the Green function Monte Carlo approach is to define a Markov process, 
yielding after a large number n of iterations a probability distribution P n (w,x) for the 
walker which determines the ground state wavefunction ipQ. To be specific in the most 
simple formulation one has: 



In the following the distribution P(w,x) is sampled by a finite number M of walkers. 
Let us first consider the simpler case M — 1. In order to define a statistical implementation 
of the matrix multiplication \x >^ H\x >, the standard approach is first to determine the 
hamiltonian matrix elements Hj x connected to x which are different from zero. Then a new 
index x' is chosen for the walker among the indices % according to the probability determined 



where bx = J2iH^ x has been introduced in order to satisfy the normalization condition 
J2iPi x = ^- This simple iteration scheme to go from a configuration a; to a new configuration 
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III. SINGLE WALKER FORMULATION 



by 



Pi,x = H i,x/ b x 



(3) 



x' is easily implemented but is not sufficient to determine stochastically the matrix-vector 
product Hx. The full matrix is a product of a stochastic matrix x and a diagonal one 

H i,x=Pi,x b x ( 4 ) 

As it is intuitively clear the diagonal matrix b x , not included in the stochastic process, is 
very easily determined by a scaling of the weight w of the walker: 

w' — > b x w (5) 

The two previous updates, the stochastic one (^) and the deterministic one (|5|) define a new 
walker w',x' in place of the "old" walker w,x, i.e. they determine a Markov process. 

At this point it is important to understand the evolution of the probability distribution 
P(w,x) after such process. A subscript n to this function P will indicate the number of it- 
erations of the Markov process. The probability evolution P n — > P n _|_ i is easily determined 
by 

P n + 1 (w',x') =J2p x > x Pn(w'/b x ,x)/b x (6) 

x 

The equation (|D allows to determine, by simple iteration, what is the probability to find 
a walker in a given configuration w, x after many steps. However the evolution P n from the 
initial distribution P is more clear and transparent in terms of its momenta over the weight 
variable w. 

Gfc n {x) = J dww k P n (w,x) (7) 

In fact it is straightforward to verify, using Eq.(||) , that : 

Gk,n + l(*') = EP x >, x bxG kjn (x) (8) 
x ' 

In particular for k = 1 the first momentum of P determines the full quantum mechanical 
information, as G\ n {x') = (H n ) i G\ q(x), implying that n (x), by Eq.(Q), converges 
to the ground state of the hamiltonian H. 

6 



By iterating several times even a single walker , the resulting configuration w,x will 
be distributed according to the ground state of H and by sampling a large number of 
independent configurations we can evaluate for instance the ground state energy: 

E = ±**> (9) 
u < w > 

where the brackets <> indicate the usual stochastic average, namely averaging over the 
independent configurations. 

This in principle concludes the GFMC scheme. However the weight w of the walker 
grows exponentially with n (simply as a result of n independent products) and can assume 
very large values, implying diverging variance in the above averages. 

In the next sections we describe in detail this problem and a way to solve it with a fixed 
number of walkers. 



IV. STATISTICAL AVERAGE DURING THE MARKOV PROCESS 

The configuration x n that are generated in the Markov process are distributed after long 
time according to the maximum right eigenstate R(x) of the matrix p x > x ( simply because 
o( x ) = Er'^Jr -t'Gq q( x ') ~ > R(x) for large n ) which, as we have seen, is in general 
different from the ground state iPq(x) we are interested in, due to the weights w n that weight 
differently the various configurations x distributed according to R(x). We are allowed to 
consider this state R(x) as the initial trial state \ipi~ > used in the power method (||), and 
that, at any Markov iteration n, the walker had weight w — 1 L iterations backward, when 
it was at equilibrium according to the distribution R(x), described before. In this way it is 
simple to compute the global weight of the walker with L power method correcting factors: 

Gn^h^n-j (10) 

Therefore, for instance , in order to compute the energy with a single Markov chain of many 
iterations, the following quantity is usually sampled: 
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(11) 



with L fixed. ]TT[ The reason to take L as small as possible is that for large L the weight 



factors G% diverge exponentially leading to uncontrolled fluctuations. In order to compute 
the variance of the factors we can simply apply what we have derived in the previous 
section and prove in few lines the exponential growth of the fluctuations of the weights. 
Using Eq.(^) it is easily found that: 



According to Eq.(§) G2 £ for large L diverges exponentially fast as \^ where A2 is the 
maximum eigenvalue of the matrix pb 2 (b is here a diagonal matrix b = 8 r >b x ), whereas 



hamiltonian matrix H = pb. It is clear therefore that we get an exponential increase of the 
fluctuations 



as in general A2 > A 2 and the equality sign holds only if the matrix b is a constant times the 
identity matrix. 

In order to overcome the problem of exponentially increasing variance, in the following 
section we will discuss a way to propagate a set of M walkers simultaneously. By evolving 
them independently, clearly no improvement is obtained for the aforementioned large fluc- 
tuations, as for this purpose it is equivalent to iterate longer a single walker. Instead, before 
the variance of the weights Wi becomes too large, it is better to redefine the set of walkers 
by dropping out the ones with a weight which is too small, and correspondingly generate 
copies of the more important ones, so that after this reconfiguration all the walkers have 
approximately the same weight. By iterating this process the weights of all the walkers are 
kept approximately equal during the simulation. This property yields a considerable reduc- 
tion of the statistical errors, as the variance of the average weight w = -k? J2i w i is reduced 
by a factor \J~M. This allows therefore a more stable propagation even for large L. 




the first momentum G^ ^ diverges as G^ ^ ~ A L , with A the maximum eigenvalue of the 



(5G L ) 2 ~ (A£ - A 2L ) 



S 



V. CARRYING MANY CONFIGURATIONS SIMULTANEOUSLY 



Given the M walkers we indicate the corresponding configurations and weights with a 
couple of vectors (w,x), with each vector component Wi, X{ i = 1, • • • , M, corresponding to 
the i th walker. It is then easy to generalize Eq. (|6]) to many independent walkers: 

P n +l(w,x) = Pn(w 1 /b Xl ,w 2 /b X2 , ■ ■ ■w M /b XM ,x' 1 ,x' 2 , ■ ■ ■ x' M ) 

x l j • h 2 i " " " ' M 

{Px 1 ,x'Px 2 , x' 2 ■ ■ -Px M , x' M ) /( b ^2 ■ ■ ■ bx M ) (12) 

If the evolution of P is done without further restriction each walker is uncorrelated from 
any other one and : 

P(wi, w 2 , ■ ■ ■ , w M , xi,x 2 ,--- x M ) = P(wi, Xi)P(w2, x 2 ) ■ ■ ■ P(w M , %m) 
Similarly to the previous case we can define the momenta over the weight variable: 

G k, n\ x ) = J dwi J dw2 "' J dwM ^ 1 ~m J Pn (—'-> 

(13) 

Since we are interested only to the first momentum of P we can define a reconfiguration 
process that change the probability distribution P n without changing its first momentum, 
and in this we follow [ |TT|j : 

P' n {wL, x_) = [ ]T G(nl, 3L] w, x)P(w, x) [dw] (14) 

X 

GU,aL;w,x) = J[ ( ^ ' ) K< ~ h ^ L ) (15) 

Hereafter the multiple integrals over all the Wj variables are conventionally shorthand by 
/ [dw\. Note that the defined Green function G is normalized / [dw'] J2 X ' G = 1. 

In practice this reconfiguration process amounts to generate a new set of M walkers 
(w'j , x'j ) in terms of the given M walkers (wj , Xj ) in the following way Each new walker w'j , x'j 
will have the same weight w = ^ 3 M J and an arbitrary configuration x'j among the possible 



old ones Xk k = 1, ■ ■ ■ , M, chosen with a probability pk = Wk/ J2j w j- It is clear that after 
this reconfiguration the new M walkers have by definition the same weights and most of the 
irrelevant walkers with small weights are dropped out. This is just the desired reconfiguration 
which plays the same stabilization effect of the conventional branching scheme. || For an 
efficient implementation of this reconfiguration scheme see the Appendix. 



VI. BIAS CONTROL 



It is well known that the control of the population size M introduces some bias in the 
simulation simply because some kind of correlation between the walkers is introduced. How- 
ever for high accuracy calculations this bias often becomes the most difficult part to control. 
In this section we can instead prove that the reconfiguration of the M walkers defined in (|HD 
does a better job. Though this reconfiguration clearly introduces some kind of correlation 
among the walkers, it can be rigorously proven that the first momentum G\ n (x) of the 
distribution of P is exactly equal to the one G'^ n (x) of P', obtained after reconfiguration. 
This means that there is no loss of information in the described reconfiguration process and 

G' hn (x) = G ln (x) (16) 

Proof: 



By definition using (|TB| ) and ([TJ 

/v . ,, ,, 

J 3 3 x, x 



G' ln (x) = J [dw] J [dw] M ' 3 G{im,d.]w,x)P{w.x) 

x,g/_\ / 

The first term in the integrand contains a sum. It is simpler to single out each term of the 
sum w' k 5 x % i /M and to integrate over all the possible variables w\ x' but w' k and x' k . It is 
then easily obtained that this contribution to G'^ n conventionally indicated as [G 1 ^ n ]k is 
given by: 



x, x k 

10 



Then by integrating simply in dw' k and summing over x' k in the previous integrand we easily 
get that [G'i n ]k = jjGi n , independent of k. Finally by summing over k we prove the 
statement fll6|) . 



VII. THE GFMC SCHEME WITH BIAS CONTROL 

Using the previous result it is easy to generalize the equations (|IT||TT|) , to many config- 
urations. It is assumed that the reconfiguration process described in the previous section is 
applied iteratively each kb steps of independent walker propagation. The index n appearing 
in the old expressions ( |10| , [TT| ) now labels the n th reconfiguration process. The measurement 
of the energy can be done after the reconfiguration when all the walkers have the same 
weight, thus in Eq. ( |T0|) : 

1 M 

£6 x n (17) 
j = l 

or, for a better statistical error, the energy can be sampled just before the reconfiguration, 
taking properly into account the weight of each single walker: 

(18) 

It is important, after each reconfiguration, to save the quantity w = jfYlfLiWj and re- 
set to one the weights of each walker. Thus the total weights G%, correspond to the 
application of Lkb power method iterations (as in Eq.||) to the equilibrium distribution 
= Gq n _ l{x) , which at equilibrium is independent of n. The value of the factor 
G^ can be easily recovered by following the evolution of the M walkers in the previous L 
reconfiguration processes and reads: 

Gi=U^n-j ( 19 ) 

3=0 

where the average weight of the walkers w has been defined previously. 

An example on how the method works for the calculation of the ground state energy of 
the HM in a 4 x 4 lattice is shown in Fig. ([!]). Remarkably our method converges very fast 
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to the exact result with few correcting factors and with smaller error bars as compared with 
the original Hetherington's scheme. Our method does not require the population control 
reconfiguration at each step, and the resulting bias is considerably reduced, in the sense 
that , without any correcting factor, our value is five times closer to the exact result, for the 
example shown in Fig. (p. 

VIII. IMPORTANCE SAMPLING 

One of the most important advantages of the Green function Monte Carlo technique 
is the possibility to reduce the variance of the energy by exploiting some information of 
the ground state wavefunction, sometimes known a priori on physical grounds. In order to 
understand how to reduce this variance, we just note that the method, as described in the 
previous sections, is not restricted to symmetric matrices, simply because we never used this 
property of the hamiltonian matrices. Following we consider not the original matrix , 
but the non symmetric one : 

where ipc is the so called guiding wavefunction , that has to be as simple as possible to 
be efficiently implemented in the calculation of the matrix elements and , as we will see, as 
close as possible to the ground state of H. 

In order to evaluate the maximum eigenvalue of H', corresponding obviously to the 
ground state of H, the coefficient b x is now given by: H' as indicated in Eq . (|Tl~|) , in this case 
b x reads: 

bx n =EMx')H x > x /M*n) (20) 
x' 

Thus if ipc is exactly equal to the ground state of H then, by definition, b Xn = Eq, inde- 
pendent of x n . This is the so called zero variance property satisfied by the method. 
Namely if the guiding wavefunction approaches an exact eigenstate of H , the method is 



12 



free of statistical fluctuations. Of course one is never in such a fortunate situation, but by 
improving the guiding wavefunction one is able to considerably decrease the error bars of 
the energy. This property, rather obvious, is very important and non trivial. 

For the application of the method to the HM we have used a Jastrow like guiding function: 
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Y,v{R-B!)S%S z R , 



\ip G >=J2s M (x)e R > R ' 



\x > 



(21) 



where \x > indicate in this case all possible spin configurations with Sr = ±4 defined on 
each site R of the I x I square lattice and with the restriction of total vanishing spin projection 
J2rSr = 0, sm(x) represents the so called Marshall sign, depending on the number N^(x) 
of spin up in one of the two sublattices sm(%) = (— l)^( x \ while the long range potential 
is given by: 

v(R) 



1- 



JqR 



1 + (cosg x + cosq y )/2 



\ 1 - (cosg x + cosg y )/2 



with \q x \ < 7T and \q y \ < ir belonging to the Brillouin zone and assuming the appropriate 
discrete values of a finite system with periodic boundary conditions. The constant 7 is the 
only variational parameter in the wavefunction, that for 7 = 1/2S* is consistent with the 



spin wave theory solution of the HM for large spin S. [12 



IX. FORWARD WALKING 

The Green function Monte Carlo can be used with success to compute also correlation 
functions on the ground state of H . In fact it is simple to compute expectation values of 
operators that are diagonal in the chosen basis, so that to a given element x of the basis 
corresponds a well defined value 0(x) = (x\0\x) of the operator. By Green Function Monte 
Carlo, as we have seen, configurations w, x distributed according to the desired wavefunction 
4>o(x), or ipQ^i/jcix) if importance sampling is implemented, are generated stochastically. 
However in order to compute < O >=< ^q|0|^q > a little further work is necessary as the 
square of the wavefunction is required to perform the quantum average. To this purpose the 
desired expectation value is written in the following form: 

13 



<jj G \H N ^OH N ' k ^ G > 



From the statistical point of view Eq.(p2|) amounts first to sample a configuration x after N' 
GFMC reconfigurations , then to measure the quantity (x\0\x) and finally to let the walker 
propagate forward for further N reconfigurations. 

In order to evaluate the stochastic average an approach similar to what done for the 



energy is clearly possible. The only change to expression (11) is to replace bxj with the av- 
erage measured quantity 0^ N = -h J2j Oj at the generation n and change the corresponding 
weight factors in (|T9|) as: 

Gn= Ef « n -j ( 23 ) 

j=-N 

where henceforth we denote with the value of the diagonal operator O on the configura- 
tion Xj of the j walker, at the iteration n. Indeed these new factors ( P3[ ) contain a further 
propagation of N reconfiguration processes as compared to the previous expression. It is 
important that both L, correcting the bias, and N, correcting the quantum average of the 
operator are finite, due to the exponential growths of the fluctuations as iV and L increase. 
On the other hand these fluctuations can be controlled by enlarging the population size M, 
and the method for M large enough remains stable. 

A further condition is however necessary in order to control the bias in the forward 
walking technique. The set of measured values Of with weight factors ( p3|) has to be modified 
after each reconfiguration process occurring in the forward direction. In practice after each 
reconfiguration it is important to bookkeep only the values of the observables that survive 
after the reconfiguration (we omit in the following the superscript n for simplicity). In other 
words after each reconfiguration 0[ = Oj(£\ f° r * = 1> ' ' ' M with the integer function j(i) 
describing the reconfiguration process in our scheme (after any reconfiguration the walker 
with index i assumes the configuration with index before the reconfiguration). 

In order to implement recursively the forward walking it is useful to store at each recon- 
figuration process the integer function j n (i) for each reconfiguration n and the values Oi of 
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the operator O for each walker. Then it is possible to compute the relevant configurations 
contributing to the operator O after N reconfiguration process by a recursive application of 
the integer functions j n , namely 0[ = O^r^^ . . - /^ . . \\\. 

It is extremely important to perform the reconfiguration process as rarely as possible since 
after each reconfiguration in the forward direction the number of different configurations 
representing an operator O decays quickly, yielding of course a larger variance. Contrary 



to the Hetherington |TT| algorithm our scheme allows the bias control without requiring the 
reconfiguration at each step. 

An example on how this scheme works is shown in Fig.(||]). As it is seen it is simple to 
reach the exact ground state average. 



X. FORMAL PROOF OF THE BIAS CONTROL IN THE FORWARD WALKING 

SCHEME 

In order to implement stochastically Eq.(|22|) we need to apply the operator O x diagonal in 
configuration space, in a stochastic sense and then follow the standard stochastic iteration 
(§) to the walker distribution P for N steps. To this purpose a walker from now on is 
identified by the triad: 

w, 7, x 

where 7 represents the actual value of the measured operator O for the walker. Its value 
can change, as we will see later on in the reconfiguration process, and in general due to the 
forward walking 7 7^ (x\0\x). Indeed only at the beginning, n = 0, of the forward walking 
iteration 7^ = Oj =< Xj|0|xj >, for i = 1, ■ ■ • M. In a probabilistic sense this is equivalent 
to consider the initial probability distribution: 

Pn=o(w, 7, X) = P (W, X) Y[ °i) 

i=l,M 

where Pq is the equilibrium distribution of the previous Markov process (||), which samples 
the ground state ipo(x). 
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With this initial condition further N forward walking steps are implemented to the 
probability distribution P, defined with the iterations in Eq.(|6]) and in Eq.(^). Then in 
order to determine the quantity ( p2| ) the following ratio is evaluated: 

< 107 > . ^ 

< O >= - — 24 

< w > 

where the brackets indicate the average over the distribution J2x I dw J d r yP(w,j,x). It is 
understood that in Eq.(J|) and Eq. ([l2]) the variables 7« remain unchanged. For instance the 
analogous of Eq. (||) will be : 

P n+1 (w',j',x') = J2p x > x Pn(w'/b x ,j',x)/b x (25) 

X 

However in order to satisfy the bias control property described in Sec. (|VTD it is necessary 
to update the 7 variables at any reconfiguration process. 

Analogously to the previous case it is easier to work with wy momenta of order k of the 
distribution P for fixed configuration x (see Eq. [7|): 



G l,n^ = J dw J 1 P n(™,7,*) (26) 

which for M 7^ 1 correspond to: 

= JdmJ i£Z { ^I2p^ Pn{m ^) (27) 

where, as usual underlined variables represent vectors whose components refer to the single 
walker index j. 

With a proof exactly analogous to the one of SeaflVT]) it is possible to show that : 

• The value of the first (u>7) momentum Gj n (x) , at the initial iteration of the forward 
walking n = 0, is equivalent to apply the operator O to the initial distribution Po(w, x), 
namely 

G l n = 0^ =°x G l,n = 0( x ) 
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the following reconfiguration process, which does not change the Markov chain of 
configurations (w,x) but modifies slightly 7, has the bias control property also for the 
wy averages: 

KfaC i_ ' = / j^G{u/_ } V, x_; w, 7, x)P(w, 7, x) [dw] \dj_ 

< = i V E ^ "^./.. / V ^ / 

The first factors in the Green function G, involving the 7's, represent the only differ- 
ence to the previous reconfiguration process flli). Thus obviously the momenta Gfc jn 
not involving the 7 variables satisfy the same bias control property of the previous 
reconfiguration process (]14|) 



As far as the (wj) momenta are concerned is possible to prove as before the mentioned 
bias control property: 

G'Ux) = G£„(s) (29) 

To this purpose, analogously to the previous case, it is convenient to single out a term 
j = k in the definition of the first wj momentum in Eq.(p7|) , and following the same route 
of Sec.flVp integrate easily the Green function over all possible variables w',7' x', but the 
variables x' k , j' k and w' k . These remaining integrations can be also performed analytically 
by first integrating in Wk, then in 7^ and finally summing over Xk- The assertion (|29|) is 
therefore proved rigorously 



XI. A "STRAIGHT" - FORWARD WALKING 

A simpler method to compute averages of general operators is obtained by Eq. ( |2~3D 
performing two independent simulations for the numerator and the denominator in the 
equation (^4|). The remarkable advantage of this technique is the possibility to measure also 
off-diagonal operators O x ^ x . After the first simulation for the evaluation of the denominator 
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we take an equilibrated walker configuration (w,x) and apply the operator O. Whenever 
the operator O is off diagonal this is not simply equivalent to scale the weight w — > wO x . 
In fact we need a stochastic approach to select only one of the configurations x' among the 
possible ones connected to x with non zero matrix element O x r jX . Whenever O x ^ x > this 
stochastic approach can be implemented with a two steps technique analogous to the one 
for the hamiltonian: 

• we first scale w by the mixed average estimate J2x' i)G( x> )Ox',x/' t l } G{ x ) 

• then we select a random new configuration x' with a probability proportional to 

^G{x')O x ^ x /i) G {x) 

Then the same reconfiguration process ([14] ) to work with a fixed number of walkers during 
the forward walking propagation can be efficiently applied. 

Finally we comment that the general operators O with arbitrary signs in the matrix 
elements can be always casted as a difference O = + — 0~ of two operators with positive 
definite matrix elements; the above method can be applied to + and 0~ separately. 

XII. DISCUSSION AND RESULTS 

In the previous sections we have described how to obtain ground state energies and 
correlation functions of some class of Hamiltonians on a finite lattice size. In this section we 
describe a successful application of this method to the HM. 

We are interested in thermodynamically converged physical quantities characterizing the 
quantum antiferromagnet, for instance the energy per site eo , the staggered magnetization 
m , the spin-wave velocity c, the spin susceptibility \. Use of a finite size scaling analysis is 
required to obtain the infinite volume limit of our data. 

We compute with our method the ground state energy per site e and the ground state 
expectation value of the spin-spin structure factor S(q) : 

S(q) = +r < H E e^ R - R '^S R ■ S R ,\% > (30) 

iVa R,R' 
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which for q equal to the antiferromagnetic momentum Q = (jr, tt) allows a finite size estimate 



of the order parameter mi = y^y • The known finite size scaling theory ||T3|,[4| in a quantum 
antiferromagnet predicts that 

a) the ground state energy per site eo has the following leading size corrections: 

e (L) = e - 1.4372 c/l 3 + ■■■ (31) 

which allows an indirect evaluation of the spin-wave velocity. 

b) Further the finite size estimate of the order parameter mi approaches its converged value 

as f 

Finally Neuberger and Ziman [|J, using a relativistic pion physics analogy, derived very 
powerful constraints on the spin-spin structure factor S(q), namely that for q — > Q it diverges 
as l/\q — Q\ with a prefactor equal to — . Using their arguments it also follows immediately 
that 

S(q) ~ X c\q\ (32) 

for \q\ — > 0. 

In spin wave theory results for the constants appearing in Eqs.(pTf) are given by: 



m = s — c' (33) 
c = 27^(1 + c /2s) (34) 
Xc = ^ (1 - c'/s) (35) 



where Cq and d can be estimated on a I x I finite size lattice: 



and 



c (/) = l-^E e ^ - 1579 ( 36 ) 

a k 



d{l) = 4r E - - - -> °- 1966 ( 37 ) 
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and are expressed in terms of the spin wave energy e k = Jl — where 7& = (cosk x + 
cos k y )/2. 

In order to improve the accuracy of the finite size scaling calculation we have systemat- 
ically compared our finite size data with the spin wave expansion on the same lattice sizes. 
]T5| This technique allows to compute explicitly the finite size corrections using the 1/S 
expansion, yielding results consistent with the previous theory for the finite size corrections. 
The advantage to use finite size spin wave theory is that it also implicitly determines all the 
subleading corrections in l/l. In this approach, the energy per site is given by: 

- Co (/)) 2 -r 



e(Z) = -2J 



(38) 



which correctly reproduces the predicted finite size scaling ( |3~1~D with the spin wave velocity 
given by Eq. (|35D and a finite size next leading contribution +2 J// 4 , appearing in the second 
order spin wave expansion. The latter term is inconsistent with a claim by Fisher, which 
probably omitted higher order contributions in his analysis. |13[ On the other hand the fit 



of the data reported in Tab,| is also in quantitative agreement with spin wave theory, also 
for this next leading contribution to the energy per site. 

By a simple Fourier transform of the finite size spin wave results for < Sr ■ So > we 
obtain S(q) = for q = 0, consistent with a singlet ground state and: 

S(Q) = N a (s - c'(Z)) 2 - l/N a + JL £ A 2 (39) 

5(g) = Ll^( a _</(/)) £ 1 ~ 7fc7 «- fc ~ ek ^ k q^Q,0 (40) 

After a simple inspection the leading behavior S(q) oc \q\ ( S(q) oc l/\q — Q\) for q — ► 
(q — > Q) is in agreement with the Neuberger and Ziman predictions, with exactly consistent 
prefactor \ c — ^75(1 — c '/ s ) ( rn 2 /x c = 2\/2s(l — d j s) ) within - g expansion. 

Now we discuss the finite size results obtained with the GFMC technique described in 
this paper. 

First we compute the ground state energy per site eo and report all the data in Tab. |. 



According to Eq. (|31|) we make a parabolic fit in l/l and obtain for the energy per site in the 
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thermodynamic limit the value reported in Tab. |I[ This value differs from previous GFMC 
estimate ||, without population control error, which is exactly removed in our method. This 
error seems to affect significantly the large size estimate of the energy as shown in the Tab.|. 
A recent paper by Sandvik |T4]| , using a completely different path integral method, reports 
energy values perfectly consistent with ours and with similar error bars. The spin-wave 
velocity is then evaluated by looking at the finite size corrections of e®. Then according to 
Eq.([31~|) we obtain: 

c /c sw = 1.12 ± 0.06 (41) 

with csw — V%J the zero order spin wave velocity. We have reached a very poor accuracy for 
this quantity, as it is determined by the subleading corrections of the energy, that in turn are 
also quite size dependent. For this quantity it is not possible to find significant differences 
from the second order spin wave result c/csw = 1.1579 (Eq.|35|), which is probably a more 
realistic and accurate value. 

The order parameter is evaluated by the forward walking technique with (Sec. [TX|) or 



without (Sec. IXT) bookkeeping the branching matrix. In the latter case the spin isotropic 



square order parameter S(Q)/N a is directly evaluated, as the method is not restricted to 
diagonal operators. As discussed in the paper the only source of systematic errors are given 
by the length of the forward walking propagation N and the number of bias correcting 
factors L for the ground state. It is understood that for N, L — > oo our method provides 
the exact finite size value for mi. The finite L error is negligible, as we have used a large 
enough number of walkers to eliminate the bias with few correcting factors. 

The most important systematic error is due to the finite N. Its drastic and controlled 
reduction as a function of iV is displayed in Fig.(^) which shows that we have achieved the 
converged values for mi within the error bars even for the largest system size N a = 16 x 16. 
Note also that the "straight" - forward method converges much more quickly as the total 
spin is conserved in this method and the convergence to the ground state is determined by 
the much larger gap in the same singlet subspace. Our data for mi are again in perfect 
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agreement with the Sandvik's ones, who was able to obtain much smaller error bars. With 
the available data we have extrapolated mi by displaying the ratio of this quantity with the 
spin- wave prediction in Fig.(^). The suppression of the finite size effects, obtained with this 
analysis is remarkably good (see table |I1|). 

Finally we have computed the structure factor with the forward walking technique 
(Sec.pC|) for several momenta. Analogously to the previous case for the order parame- 
ter, we have studied the ratio of the QMC data with the spin wave prediction given by 
Eq.fl39|). This calculation shows that the spin wave expression (|39|) is particularly accu- 
rate close to q ~ Q but there is some deviation for small momenta. From Eq . (|39|) the 
small q limit of the structure factor can be computed analytically in spin wave expansion: 
S(q) = \q\/2V2(s - d) = 0.1073|g|. The ratio between our QMC data and the second order 
spin wave prediction is shown in Fig. ([5]) and for small q it approaches the value 1.045 ±0.01, 
rather independently from the system size. This yields, by Eq. (|32|), a direct determination of 
the product of the spin wave velocity and the susceptibility: \ c = 0.1122±0.001. This value 
is very much in disagreement with Sandvik's one (JTJJ who predicted x c = 0.1046 ± 0.002, 
about four error bars out from our direct and more accurate result. This discrepancy is prob- 
ably due to a very difficult infinite size extrapolation of x [PI and c, whereas the product 



X x c calculated by means of S(q) for small momenta, seems rather well behaved, as shown in 
Fig.([|). In the same figure the prefactor around q ~ Q is in quite reasonable agreement with 
the expected one (full big dot) m 2 /xc obtained with our independent measures of x c (slope 
S q for q — > 0) and m 2 (S(q)/N a for q = Q), considering also that this function, according to 
the Neuberger and Ziman's theory, should jump discontinuously at q = Q, where it assumes 
a value determined only by the the order parameter m 2 . The shape of the curve displayed 
in Fig,(|5|) is clearly consistent with the predicted singularity. 

Our results therefore represent a direct confirmation about the internal consistency of 
the Neuberger and Ziman's theory. 

The values for the physical quantities extrapolated in the thermodynamic limit is sum- 
marized in table (|Hl|). 
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XIII. CONCLUSIONS 



We have described in detail a straightforward implementation of the Green Function 
Monte Carlo scheme on a lattice hamiltonian without need of the standard branching process. 
The extension of such a scheme to continuum systems such as He 4 is straightforward. Indeed 
the present algorithm works at fixed number of walkers and we have shown that all sources of 
systematic error can be controlled with a rigorous approach both for computing the ground 
state energy and for computing ground state correlation functions. The possibility to work 
with a limited number of walkers is extremely important for high accuracy calculations. 
This in fact requires quite long simulations to decrease the statistical errors and, with the 
standard approach f2|, one easily exceeds the maximum number of walkers for the available 
computer memory. 

Our reconfiguration scheme is not restricted to work at fixed number of walkers. In fact 
our proofs in Sections |V]],[X] can be readily generalized when the number M' of outgoing 
walkers x'j (with unit weight) is different from the number M of the incoming ones Wj Xj. 
Thus a standard branching scheme between two consecutive reconfigurations can be also 
applied, and the method can be used only each time the population of walkers reaches an 

M 

exceedingly large or small size. At each stochastic reconfiguration the simple factor J2 Wj/M' 
corrects exactly the bias of the described size population control. This maybe a more efficient 
implementation of our method, with, however, a rather more involved algorithm. 

This scheme represents a more practical implementation of the Hetherington idea to 
work at fixed number of walkers as: 

i) it is not required to apply the reconfiguration process at each Markov iteration. Indeed 

our scheme coincides with the Hetherington one in this limit. 

ii) the same idea to control the bias at fixed number of walkers was extended to the "forward 

walking technique" which is important to calculate efficiently correlation functions on 
the ground state. 
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iii) contrary to the conventional believe we have shown that there are no basic difficulties to 
compute "off diagonal" correlation functions with Green Function Monte Carlo, using 
a simple forward walking technique (SecjXip, which turned out to be very efficient for 
determining the order parameter in the HM. 

We have applied these methods to compute very accurately the ground state energy per 
site of the HM, the spin-spin structure factor and the antiferromagnetic order parameter, 
whose infinite size values are shown in the table ([LLJ|) . 

We have obtained very good agreement with finite size spin wave theory, which allows 
a very well controlled finite size scaling (see Fig|J). We believe that the reported accuracy 
gives a very robust confirmation about the existence of AF long range order in the 2D HM. 

We finally show in Fig.(^) a comparison of our QMC prediction for S(q) with the available 
experimental data on the stechiometric La2Cu04 Mott insulator for momenta close to the 
AF wavevector: q ~ (tt, 7r). The agreement is remarkably good considering also that 

i) the experiments have been performed at a temperature above the Neel temperature, so 

that no true long range order exists, despite the very long correlation length. As a 
consequence, for a good comparison between experiments and theory, one should take 
into account the smearing of the 5 function contribution at Q = (tt, 7r), to be added 
to the theorethical ground state prediction. 

ii) there are no fitting parameters in the present comparison 

Therefore the Copper Oxygen planes of La2Cu04, planes which become High Tempera- 
ture superconductors upon finite hole doping, are well described by the nearest neighbor 
Heisenberg model. 
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APPENDIX: PRACTICAL RECONFIGURATION PROCESS 

In this appendix we follow the notations of Sec.flVTD to describe an efficient implementa- 
tion of the reconfiguration process needed to stabilize the proposed GFMC method. 

The new walkers x\ after reconfiguration are chosen among the old ones x k , with proba- 
bility pk- We divide the interval (0, 1) in M subintervals with length pk from the leftmost to 
the rightmost as k increases from 1 to M. Then we generate M pseudo-random numbers Zj 
for i — 1, • • • M and sort them, baring in mind that the index i labels the walker x\ after the 
reconfiguration. We save therefore the random permutation i(k) k = 1, . . . M corresponding 
to the described sorting, permutation that determines > z^ k y An efficient sorting 

algorithm takes the order of M log M operations, thus is not time consuming. 

The next step p is to make a loop over the sorted index k, giving monotonically increasing 
Zj(fc), and to select as new configuration x'^ the one Xj, among the old configurations, such 
that ZiOA belongs to the interval of length pj. Note that the index function % = 1, • • • M 
contains all the information required for the forward walking technique described in Sec. (|LX| ) . 

After the described process some of the old configurations may appear in many copies, 
while others disappear. This happens also if the distribution pj is uniform pj ~ 1/M, yielding 
clearly some loss of information in the statistical process. A better way to implement 
the reconfiguration, without loosing information and without introducing any source of 
systematic error, is obtained by the following simple change, After the generation of the 
random permutation i(k) a new set of numbers uniformly distributed in the interval (0, 1) 
is defined : 

zt = (Z+(i-l))/M 

for % — 1, • • • M, where £ is another pseudorandom number in (0, 1). This set of numbers z~i, 
now uniformly distributed in the interval (0, 1), is then used to select the new configurations, 
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yielding a more efficient implementation of the described reconfiguration process. 
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TABLES 

TABLE I. Energy per site of the HM in the square lattice I x I . In this work the infinite size 
extrapolation is obtained with a parabolic fit Eq(1) = Eq{oo) + a/Z 3 + b/l 4 for all size with I > 8. 
( a = —2.275 ± 0.12 and , b = 1.64 ± 0.9). The numbers in parenthesis represent error bars in the 
last digits. 
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Eo i 


E [14] 


-Eo This work 


6 


-0.678871(8) 


-0.678873(4) 


-0.6788721(28) 


8 


-0.673486 (14) 


-0.673487(4) 


-0.673483(8) 


10 


-0.671492(27) 


- 0.671549(4) 


-0.671554(6) 


12 


-0.670581(49) 


-0.670685(5) 


-0.670678(5) 


14 




-0.670222(7) 


-0.670223(8) 


16 


-0.669872(28) 


-0.669976(7) 


-0.669977(8) 


oo 


-0.66934(3) 


-0.669437(5) 


-0.669442(26) 



TABLE II. Staggered magnetization of the HM in the square lattice I x I with side I computed 
with the forward walking technique. The numbers in parenthesis represent error bars in the last 
digits. The extrapolated values of our GFMC data for N a — > oo are obtained by the fit shown in 
Fig! 
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mi § 


mi [14] 


mi This work, sec. IX 


This work, sec. XI 


6 


0.4581(2) 


0.458074(3) 


0.4583(3) 


0.4579(1) 


8 


0.420(1) 


0.421709(9) 


0.4212(6) 


0.4217(1) 


10 


0.397(3) 


0.399214(9) 


0.3988(8) 


0.3991(6) 


12 


0.378(14) 


0.38400(1) 


0.380(2) 


0.3834(6) 


16 




0.3647(1) 


0.361(9) 


0.3642(6) 


oo 


0.3075(25) 


0.3070(3) 


0.3058(12) 


0.3077(4) 
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TABLE III. Infinite size estimates of the various ground state quantities of the HM discussed 



in the paper. The value of m is obtained by fitting the more accurate data of [14 as shown in 
FigJ^. The value of x c is computed from the data of Fig.([|), corresponding to the largest size and 
smallest momentum. 



eo 



-0.669442(26) 



m 



0.3075 ± .0002 



Xc 



0.1122 ± 0.001 



28 



FIGURES 

FIG. 1. Energy per site in a 4 x 4 Heisenberg cluster. The continuous line is the exact result 
and the dotted (dashed) line connects GFMC data as a function of the number L of correcting 
factors within our scheme (Hetherington's one) described at the end of the appendix. The number 
of walkers in this case was M = 10 and the reconfiguration scheme was applied each four (k p = 4 
in the text) iterations , while in the Hetherington's scheme k p = 1. The guiding wavefunction is 
given in Eq.(^) with 7 = 1.2, and all the data are obtained with the same amount of computer 
time. The insect is an expansion of the bottom-left part of the picture. 

FIG. 2. Plot of the squared antiferromagnetic order parameter mf, in a 4 x 4 Heisenberg cluster 
as a function of the number N of forward walking reconfigurations. The continuous line indicates 
the exact result. The number of walkers was fixed to M = 20 with k p = 5. The guiding function 
was the one referenced in the Fig.Q. With the present technique the z component of mi can be 
measured on each sampled configuration ; mf = j^J2r(~^) R Sr an d spin isotropy is used to 
determine mf = 3 < (mf) 2 >. 

FIG. 3. Staggered magnetization for increasing lattice sizes (from top curve to bottom curve) 



as a function of the forward walking iteration number N computed with the method of section IX 
(a) and with the method of section [X| (b). The number of walkers for each lattice size and figure 
(a) is M = 1000, 2000, 3000, 3000, 3000 with k p = 30, 50, 60, 80, 80 and I = 6, 8, 10, 12, 16 
respectively, while the guiding wavefunction is given by Eq.([2l]) with 7 = 1.125. In picture (b) the 
number of walkers is M = 1000, 2000, 2000, 2000, 2000 and kp = 10, 20, 20, 20, 25 

FIG. 4. Plot of the ratio between staggered magnetization, computed by forward walking 
GFMC, and the spin- wave staggered magnetization as a function of 1/7. The circles refer to 
the "straight" forward walking technique ( Sec. [X|), and the squares to the data reported in Ref. 
]j~4|| . The lines are weighted linear least square fit of the data. The extrapolations for I — ► 00 are 
displayed in Tab. || and in Tab. II I 
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FIG. 5. Plot of the ratio between the spin-spin structure factor S(q), computed 
by forward walking GFMC, and the second order spin wave estimate Sswil) ( see 
Eqj4C|) . The Structure factor has been evaluated for L = 6,8,10,12 over the path 
r = (0,0) M = (vr,0) -> Y = (vr,7r) — > T = (0,0). The big dot at the Y point is the 

S(q) 

Ssw(q) ' 



expected q — > Y limit of c S( ~ 9 } n \ , by using the values for m 2 and x c reported in Tab. ID 



FIG. 6. Comparison between QMC data (full squares) and experimental data (full dots) for 
the static magnetic structure factor. The continuous line connecting the QMC data is the second 
order spin wave prediction (see Eq. ^), almost exact in this scale, whereas the line connecting the 
experimental data represents the corresponding fit reported in M. 
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